knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
fig.height = 4, fig.width = 7)
library(terra)
library(sf)
library(data.table)
library(tidyverse)
library(here)
library(oharac)
source(here('common_fxns.R'))
Using functional entities defined in a prior script based on a suite of categorical traits, in combination with spatial distribution of species from AquaMaps and IUCN, calculate functional diversity metrics based on Mouillot et al (2014) and code from Sebastien Villeger: functional richness, functional vulnerability, functional redundancy, and functional over-redundancy.
Included species are:
Here we use the IUCN species range maps rasterized to 10 km Mollweide, and AquaMaps species distributions reprojected to 10 km Mollweide. Priority: IUCN species maps, then AquaMaps with occur_cells ≥ 10. These map filenames are already set up in the assemble_spp_info_df() function.
Read in functional entity assignments and join to spatial ranges. Per cell per FE, calculate the number of species for later calculation into functional richness, redundancy, vulnerability. Break the process into chunks to avoid overloading memory.
process_spp_to_fe <- function(fe, spp_maps) { ### fe <- fe_vec[1]
dt_out <- spp_maps %>%
filter(fe_id == fe) %>%
data.table() %>%
.[ , .(n_spp_fe = length(unique(species))),
by = .(cell_id, fe_id)]
return(dt_out)
}
fe_summary_file <- here_anx('func_entities', 'fe_cell_summary.csv')
# unlink(fe_summary_file)
if(!file.exists(fe_summary_file)) {
spp_info_df <- assemble_spp_info_df(fe_only = TRUE, vuln_only = TRUE)
### table(spp_info_df %>% select(species, len) %>% distinct() %>% .$len)
# huge lg med ml sm tiny vlg vsm
# 406 1370 3595 2555 4179 2174 906 5959
### Chunk out and save smaller files to avoid loading in EVERY spp map
### all at once!
n_chunks <- 10
spp_to_chunks <- spp_info_df %>%
select(species) %>%
distinct() %>%
mutate(chunk = rep(1:n_chunks, length.out = n()))
spp_clean <- spp_info_df %>%
select(species, fe_id, map_f, taxon) %>%
distinct() %>%
left_join(spp_to_chunks, by = 'species')
spp_fes <- spp_clean %>%
select(species, fe_id) %>%
distinct()
tmp_chunk_filestem <- here('tmp/fe_summary_chunk_%s.csv')
for(i in 1:n_chunks) { ### i <- 1
tmp_chunk_f <- sprintf(tmp_chunk_filestem, i)
if(file.exists(tmp_chunk_f)) {
message('Temp file ', basename(tmp_chunk_f), ' exists... skipping!')
next()
}
message('Assembling species maps and binding FEs for chunk ', i, ' of ', n_chunks, '...')
spp_chunk <- spp_clean %>% filter(chunk == i)
chunk_spp_maps <- collect_spp_rangemaps(spp_vec = spp_chunk$species,
file_vec = spp_chunk$map_f,
idcol = 'species',
) %>%
dt_join(spp_fes, by = 'species', type = 'left') %>%
distinct()
fe_vec <- chunk_spp_maps$fe_id %>% unique() %>% sort()
message('... summarizing FR map for chunk ', i, ': ', length(fe_vec), ' FEs in ',
nrow(chunk_spp_maps), ' cell observations...')
chunk_fe_sumlist <- parallel::mclapply(fe_vec, mc.cores = 40,
FUN = process_spp_to_fe,
spp_maps = chunk_spp_maps)
if(check_tryerror(chunk_fe_sumlist)) {
stop('Try-errors detected in FE summary loop!')
}
message('... binding chunk and writing chunk summary to ', basename(tmp_chunk_f), '...')
chunk_fe_summary <- chunk_fe_sumlist %>% data.table::rbindlist()
fwrite(chunk_fe_summary, tmp_chunk_f)
}
message('... reading temp chunks, binding, and summarizing by fe and cell...')
fe_tmp_fs <- list.files(dirname(tmp_chunk_filestem),
pattern = 'fe_summary_chunk', full.names = TRUE)
fe_summary_list <- parallel::mclapply(fe_tmp_fs, FUN = fread, mc.cores = 10)
fe_summary <- rbindlist(fe_summary_list) %>%
.[ , .(n_spp_fe = sum(n_spp_fe)),
by = .(cell_id, fe_id)]
message('... writing summary to ', basename(fe_summary_file), '...')
fwrite(fe_summary, fe_summary_file)
### unlink(fe_tmp_fs)
}
Calculate functional diversity metrics from IUCN and AquaMaps summaries of functional entity membership per cell.
This set of calculations replicates the results of Villeger’s function but without resorting to a presence/absence matrix, instead just relying on the original AquaMaps species-cell table.
fe_metrics_file <- here_anx('func_entities/fe_metrics_per_cell.csv')
### unlink(fe_metrics_file)
if(!file.exists(fe_metrics_file)) {
# system.time({
fe_sum_total <- data.table::fread(fe_summary_file)
### there are 6.5 million cells. Chop into 500k instances and mclapply to summarize
# cell_id_vec <- 1:ncell(raster::raster(here('_spatial/ocean_area_mol.tif')))
chunk_size <- 500000
n_chunks <- 6.5e6 / chunk_size
message('Calculating species richness map...')
nspp_list <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
FUN = function(n) { ### n <- 1
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
nspp_sum <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
data.table() %>%
.[ , .(n_spp = sum(n_spp_fe, na.rm = TRUE),
n_fe = n_distinct(fe_id)),
by = 'cell_id']
})
if(check_tryerror(nspp_list)) stop('Try errors detected!')
nspp_df <- nspp_list %>%
data.table::rbindlist()
message('Calculating functional vulnerability map...')
fv_list <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
FUN = function(n) { ### n <- 1
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
f_vuln <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
oharac::dt_join(nspp_df, by = 'cell_id', type = 'left') %>%
data.table() %>%
.[ , .(f_vuln = sum(n_spp_fe == 1) / n_fe,
f_wvuln = mean(0.5^(n_spp_fe - 1)),
f_red = mean(n_spp_fe)),
by = .(cell_id, n_fe, n_spp)]
})
if(check_tryerror(fv_list)) stop('Try errors detected!')
f_vuln_df <- fv_list %>%
data.table::rbindlist()
message('Calculating functional overredundancy map...')
f_overred_list <- parallel::mclapply(1:n_chunks, mc.cores = n_chunks,
FUN = function(n) { ### n <- 1
cell_id_min <- (n - 1) * chunk_size + 1
cell_id_max <- n * chunk_size
f_overred <- fe_sum_total %>%
filter(between(cell_id, cell_id_min, cell_id_max)) %>%
oharac::dt_join(f_vuln_df, by = 'cell_id', type = 'left') %>%
data.table() %>%
.[ , f_over := ifelse(n_spp_fe > f_red, n_spp_fe - f_red, 0)] %>%
.[ , .(f_overred = sum(f_over) / sum(n_spp_fe)),
by = 'cell_id']
})
if(check_tryerror(f_overred_list)) stop('Try errors detected!')
f_overred_df <- f_overred_list %>%
data.table::rbindlist()
message('Combining functional diversity metrics maps...')
fe_metrics_df <- f_vuln_df %>%
oharac::dt_join(f_overred_df, by = 'cell_id', type = 'full')
write_csv(fe_metrics_df, fe_metrics_file)
}
Create and save rasters of functional diversity metrics
fe_metrics_df <- data.table::fread(fe_metrics_file)
### note that while there are 3,684,273 cells in this dataframe,
### 4273 of those are the Caspian Sea, so when those are dropped,
### it results in 3,680,000 exactly... ugh, weird
n_spp_rast <- map_to_mol(fe_metrics_df, which = 'n_spp')
n_fe_rast <- map_to_mol(fe_metrics_df, which = 'n_fe')
f_vuln_rast <- map_to_mol(fe_metrics_df, which = 'f_vuln')
f_wvuln_rast <- map_to_mol(fe_metrics_df, which = 'f_wvuln')
f_red_rast <- map_to_mol(fe_metrics_df, which = 'f_red')
f_overred_rast <- map_to_mol(fe_metrics_df, which = 'f_overred')
writeRaster(n_spp_rast, here('_output/func_entities/n_spp.tif'),
overwrite = TRUE)
writeRaster(n_fe_rast, here('_output/func_entities/n_fe.tif'),
overwrite = TRUE)
writeRaster(f_vuln_rast, here('_output/func_entities/f_vuln.tif'),
overwrite = TRUE)
writeRaster(f_wvuln_rast, here('_output/func_entities/f_wvuln.tif'),
overwrite = TRUE)
writeRaster(f_red_rast, here('_output/func_entities/f_red.tif'),
overwrite = TRUE)
writeRaster(f_overred_rast, here('_output/func_entities/f_overred.tif'),
overwrite = TRUE)
Examine distributions of each raster.
Plot each raster. NOTE: These no longer include spp that are missing from vulnerability calculations - so the maps in _output/nspp_maps (vuln \(\cap\) FE) should be identical to these in output/func_entities.